library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.2
library(corrplot)
## corrplot 0.84 loaded
library(qvalue)
## Warning: package 'qvalue' was built under R version 3.5.2
IBD_ldsc = get(load("/Users/kushaldey/Documents/singlecellLDSC/output/Healthy_gene_score_Feb1_plus_All_IBD.rda"))
dim(IBD_ldsc)
## [1] 51 64  3  4

tau-star

Roadmap Enhancers linked to Genes

tau_table = IBD_ldsc[,,3,1]
ptau_table = IBD_ldsc[,,3,2]
qtau_table = matrix(qvalue(as.vector(ptau_table))$qvalues, nrow = nrow(tau_table), ncol = ncol(tau_table))
ptau_table2 = matrix(qvalue(as.vector(ptau_table))$pvalues, nrow = nrow(tau_table), ncol = ncol(tau_table))
tau_table[which(qtau_table > 0.1)] = 0
#tau_table[tau_table < 0] = 0
corrplot(tau_table, is.corr = F, tl.cex = 0.7, tl.srt = 45)

apply(tau_table, 2, mean)[order(apply(tau_table, 2, mean), decreasing = T)][1:10]
##                      PASS_IBD_deLange2017 
##                                 0.3114454 
##          UKB_460K.blood_RBC_DISTRIB_WIDTH 
##                                 0.2301954 
##  UKB_460K.disease_HYPOTHYROIDISM_SELF_REP 
##                                 0.1990420 
##             UKB_460K.blood_PLATELET_COUNT 
##                                 0.1922581 
## UKB_460K.disease_ALLERGY_ECZEMA_DIAGNOSED 
##                                 0.1760991 
##                 PASS_Rheumatoid_Arthritis 
##                                 0.1670322 
##                PASS_Alzheimers_Jansen2019 
##                                 0.1575800 
##                UKB_460K.blood_WHITE_COUNT 
##                                 0.1480575 
##                  UKB_460K.blood_RED_COUNT 
##                                 0.1476459 
##        UKB_460K.biochemistry_TotalProtein 
##                                 0.1389440

100kb linked to Genes

tau_table = IBD_ldsc[,,3,1]
ptau_table = IBD_ldsc[,,3,2]
qtau_table = matrix(qvalue(as.vector(ptau_table))$qvalues, nrow = nrow(tau_table), ncol = ncol(tau_table))
ptau_table2 = matrix(qvalue(as.vector(ptau_table))$pvalues, nrow = nrow(tau_table), ncol = ncol(tau_table))
tau_table[which(qtau_table > 0.1)] = 0
corrplot(tau_table, is.corr = F, tl.cex = 0.7, tl.srt = 45)

Enrichment

Roadmap Enhancer

E_table = IBD_ldsc[,,3,3]
pE_table = IBD_ldsc[,,3,4]
qE_table = matrix(qvalue(as.vector(pE_table))$qvalues, nrow = nrow(E_table), ncol = ncol(E_table))
pE_table2 = matrix(qvalue(as.vector(pE_table))$pvalues, nrow = nrow(E_table), ncol = ncol(E_table))
E_table[which(qE_table > 0.01 & pE_table2 > 0.001)] = 1 
E_table[E_table < 1.1] = 1
Edif_table = E_table - mean(E_table)
Edif_table[Edif_table < 0] = 0
corrplot(Edif_table, is.corr = F, tl.cex = 0.7, tl.srt = 45)

apply(E_table, 2, mean)[order(apply(E_table, 2, mean), decreasing = T)][1:10]
##                     PASS_IBD_deLange2017 
##                                12.133796 
##         UKB_460K.blood_RBC_DISTRIB_WIDTH 
##                                10.034983 
##       UKB_460K.biochemistry_TotalProtein 
##                                10.019738 
## UKB_460K.disease_HYPOTHYROIDISM_SELF_REP 
##                                 9.557804 
##                                 PASS_LDL 
##                                 9.482519 
##            UKB_460K.blood_PLATELET_COUNT 
##                                 9.385778 
##               PASS_Alzheimers_Jansen2019 
##                                 9.143306 
##                 UKB_460K.blood_RED_COUNT 
##                                 8.696251 
##                PASS_Rheumatoid_Arthritis 
##                                 8.674550 
##               UKB_460K.blood_WHITE_COUNT 
##                                 8.457306

100kb

E_table = IBD_ldsc[,,1,3]
pE_table = IBD_ldsc[,,1,4]
qE_table = matrix(qvalue(as.vector(pE_table))$qvalues, nrow = nrow(E_table), ncol = ncol(E_table))
pE_table2 = matrix(qvalue(as.vector(pE_table))$pvalues, nrow = nrow(E_table), ncol = ncol(E_table))
E_table[which(qE_table > 0.01 & pE_table2 > 0.001)] = 1 
E_table[E_table < 1.1] = 1
Edif_table = E_table - mean(E_table)
Edif_table[Edif_table < 0] = 0
corrplot(Edif_table, is.corr = F, tl.cex = 0.7, tl.srt = 45)

apply(Edif_table, 2, mean)[order(apply(E_table, 2, mean), decreasing = T)][1:10]
##                UKB_460K.blood_RBC_DISTRIB_WIDTH 
##                                       0.6817303 
##                   UKB_460K.blood_PLATELET_COUNT 
##                                       0.6014767 
##                      PASS_Alzheimers_Jansen2019 
##                                       0.5986725 
##                                        PASS_LDL 
##                                       0.6014308 
##               UKB_460K.biochemistry_Cholesterol 
##                                       0.5424124 
##              UKB_460K.biochemistry_TotalProtein 
##                                       0.5272918 
##       UKB_460K.biochemistry_AlkalinePhosphatase 
##                                       0.5150232 
##            UKB_460K.biochemistry_TotalBilirubin 
##                                       0.4795640 
##                            PASS_IBD_deLange2017 
##                                       0.4669224 
## UKB_460K.biochemistry_AspartateAminotransferase 
##                                       0.4445825